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I. INTRODUCTION 



Molecular dynamics (MD) simulations of solute-solvent systems in chemistry and biology 
require accurate computation of electrostatic forces in order to obtain meaningful results. For 
practical purposes, computational efficiency is also essential, and various formulations exist that 
strive to achieve a balance between these two requirements. The explicit solvent methods simulate 
behaviour of each single solvent molecule which may be prohibitively expensive for a system of 
reasonable dimensions. In addition to having high computational costs, explicit solvent methods 
are usually tailored for reproducing one of the many physical properties of the solvent and 
therefore may not be well suited for a general description of solute-solvent systems (see [l|, 0, 0] 
for reviews and performance analyses). 

The alternative approach is to treat the solvent as a dielectric continuum, and the solute as a 
different dielectric object in the solvent. The dielectric properties of the solvent and the solute 
usually serve as parameters of the model. In the literature this scheme is known as the implicit or 
continuum solvent method (for reviews see (4,0, f], 0, 0]). Computations based on these methods 
are inherently faster while comparable in accuracy with those using explicit methods, at least 
in the situations when interactions between solute and solvent molecules can be neglected. For 
reasons of computational efficiency, many of the implemented implicit solvent methods make use 
of assumptions which prevent improvement in accuracy even as computational resources increase. 
The so-called generalized Born model is a good example of such uncontrolled approximations (see 
{9] for a discussion). 

To achieve controllable accuracy, we have recently proposed a novel scheme Q based on 
determining surface charges satisfying the displacement field boundary condition. With this 
scheme, one can achieve any level of accuracy permitted by the available computing power, while 
remaining computationally more efficient than explicit solvent methods. The main idea is to 
treat the induced surface charges at the boundaries as the variables to be solved for. This makes 
the potential, expressed directly in terms of the induced surface charge density, continuous at 
the boundary. Therefore only the displacement field boundary condition remains, and it leads 
to a set of algebraic equations for the surface charge densities. The potential is obtained at no 
additional cost. 

One of the seeming oversimplifications in the implicit solvent methods is the assumption of a 
sharp boundary between the solute and solvent. It is known, for example, that the solute (e.g., 
proteins) may stro ngly interact with the surrounding solvent molecules producing the so-called 
hydration layer (s) 10] . To determine electrostatic forces acting on a protein coated with such 
hydration layers, one needs to find induced charges in a spatially varying dielectric medium. 
In this paper we develop a rigorous framework, based on functional minimization, for handling 
spatially varying dielectrics. 

Functional variation is a powerful approach in modern physics. Despite common use in quan- 
tum electrodynamics, variational techniques in classical electrostatics are relatively rare and 
focus mainly on boundary value problems for linear dielectrics ll|. It has long been a textbook 
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fact that the true electrostatic potential minimizes the system's energy for a given configura- 
tion of charges 121 ] . A suitable energy functional can be constructed in general for any system 
of continuous media including systems with inhomogeneous and nonlinear dielectric properties. 
For instance, free energy functionals became an important tool in description of electrolyte solu- 
tions within the mean-field (Poisson-Boltzmann) approach (see a recent paper 13j and references 
therein). 

From our viewpoint the electrostatic potential is not the best choice for a minimization variable 
as it contains information about both the cause and effect, i.e. the source and induced charge 
densities. Moreover constitutive relations must be assumed (as in 14(). And finally this approach 
depends on prior knowledge of the Green's function with boundary conditions suitable for the 
given problem. In contrast, we use the polarization as the fundamental function as was proposed 
by Marcus over fifty years ago 15[, albeit with a different functional. The constitutive relations 
are then obtained as a result of minimization of the energy functional. The only boundary 
condition needed is that the potential goes to zero sufficiently rapidly (like inverse of the distance) 
at large distances. 

In Marcus's formulation 15|, the electric field and electric polarization were strongly motivated 
as the vectors defining the electric state of the system. This formulation was aimed at the 
processes (charge transfer chemical reactions) which happen on a much shorter time scale than 
the molecular rearrangement in response to the changing electric field. Marcus attempted to 
deal with this problem by dividing the polarization into a fast reacting part that is proportional 
to the local electric field and slowly reacting part that is not a function of the local electric 
field. As a result, the free energy functional derived by Marcus contains several electric fields 
and polarizations of various origins. 1 

A physically sound free energy functional was proposed by Felderhof 16(] in the context of a 
discussion of thermal fluctuations of the polarization and magnetization in dielectric magnetic 
media. However, a free energy of this type seems not to have been adopted for calculation of 
electrostatic fields until recently. An example of numerical implementation using Felderhof's 



scheme can be found in 17J. There, the polarization vector field was expanded in a plane 



waves basis set. The energy functional is then an ordinary function of the expansion coefficients 
which, in turn, become the variational parameters of a standard multidimensional optimization 
problem. Fast Fourier transforms were used to go from the real space to the reciprocal space 
representations. 

An approach close to Felderhof's scheme was also taken in 18[ where a thermodynamic func- 



If there is a true separation of time scales between various portions of the electrical response, these excess fields 
should be eliminated by a proper classification of the charges in the system: charges that respond rapidly and 
whose redistribution is a function of the local electric field contribute to the polarization, while charges that 
respond slowly are part of the so-called free charge distribution. However, if one were allowed to combine the 
induced charge due to fast-responding polarization with the frozen free charges, Marcus's functional becomes 
identical to ours. 
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tional was constructed with the polarization as the independent function. However, the tech- 
niques used there are suitable for systems with sharp boundaries only (the susceptibility is not 
considered to be a function of coordinates, but is rather treated as a piecewise constant). 

In this paper we construct an energy minimization scheme suitable for a rigorous treatment 
of systems with spatially varying dielectric functions, be they linear or nonlinear. In the case 
of linear dielectrics, our functional is equivalent to that proposed by Felderhof [16j ]. In section 
UTI we give the details of the formulation and describe a systematic protocol for obtaining the 
total charge density. To show the versatility of the scheme we apply it in section II I II to systems 
with sharp boundaries for which the exact solutions (or the exact equations governing the exact 
solutions) are known. In section IIVI we present numerical results for the case of two interacting 
dielectric charged spheres (solutes) placed in a dielectric solvent. We discuss the differences in 
force and energy between the situations with sharp and smooth boundaries. Finally we conclude 
with a discussion assessing the usefulness of the method. Electrostatic CGS units are used 
throughout. 



II. FUNDAMENTAL FORMULATION 

Polarization is the response of a dielectric medium to an applied electric field. The phe- 
nomenon is usually visualized as the appearance of an induced dipole moment due to a small 
shift in the relative positions of the positive and negative charge centers at the atomic scale 
[igl | . The shift may be either translational or rotational or both, depending on the quantum 
mechanical and electromagnetic interactions at the atomic level. The applied electric fields must 
be weak enough not to split the atoms or molecules into their constituents. The system is in a 
state of equilibrium under the external electromagnetic and the intrinsic restoring forces. 

Quantitatively, polarization P(r) is the density of induced dipole moment at location r. This 
density in classical electrodynamics is defined through averaging of dipole moments of con- 
stituent atoms/molecules in a small volume centered around r. The amount of polarization 
depends on the applied force and the susceptibility of the medium to such forces. Determination 
of the susceptibility of the medium (or rather the intrinsic restoring force in the medium) is the 
subject of quantum mechanics rather than classical electrodynamics. Polarization is thus a clas- 
sical/macroscopic variable summarizing quantum mechanical effects at the atomic/microscopic 
level. Therefore, we choose the polarization vector field P(r) and electric field E(r), in contrast 
to the more commonly used pair E(r) and D(r), as our fundamental variables. This choice 
provides a simpler connection to the parameters determined in microscopic physics. 

We express the energy as a functional U[P] 

U[P] = U C [P] + W[P], (1) 

where E/ic[P] is the electrostatic energy of interaction of all charges present in the system, and 
W[P] is the energy required to create the given polarization vector field P(r). 
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From simple considerations it can be shown 19|, |20(] that the variation of polarization in the 
vicinity of a point is equivalent to the presence of an induced charge density p;(r) = —V ■ P(r). 
Therefore, the total charge density Pt( r ) m the medium is a sum of the free charge density Pf(r) 
and Pi(r): 

pt(r) = p f (r)+ Pi (r). (2) 

Then 2 



[Pf(r) 



V ■ P(r) 



1 



[pf(r') - V ■ P(r')] drdr'. 



(3) 



Note first that we do not include any separate term for induced surface charges as was done in 



some of the earlier formulations of functional minimization |15|, [18] . The volume charge density is 
the most general form of charge density possible. Secondly, ([3]) is the Coulomb energy in vacuum 
and hence quite fundamental as opposed to the form with the dielectric constant of the material 
in the denominator used in some of the earlier works [3] . 

The work functional W[P] should contain the intrinsic self interaction of the polarization 
vector field. Here we consider only local contact terms for the intrinsic interactions. Noting that 
the energy functional is a scalar and assuming P <-> — P symmetry, one can write the general 
work functional W[P] as a polynomial expansion in even powers of P (or the components Pi). 
Thus we may write, 



W\P] 



Pi 



1 



x( r ) 



Pj + PPj 



Mr) 



PkPi 



ijkl 



dr, 



(4) 



where the interaction tensors etc. describe the linear and nonlinear dielectric properties 

of the media, isotropic or anisotropic (summation over repeated indices is assumed). The effective 
dielectric properties of the medium at the macroscopic level are now contained in these quantities. 

We emphasize that UJP] is the actual energy functional unlike various other functionals pro- 
posed in the literature 12, 3, 0, 0] which yield the energy or free energy of the system only at 
equilibrium. The equilibrium distribution of polarization (as well as induced charge distribution) 
can be obtained by minimizing this energy functional with respect to the polarization. For any 
given external charge distribution and spatially varying dielectric susceptibilities one can obtain 
the solution analytically or numerically. 

We may truncate the series in (j3J) at an order suitable for the problem at hand. For example, 
if the field is very weak we can retain only the quadratic term which corresponds to the case of 
linear dielectrics (isotropy is also assumed for the sake of simplicity of presentation): 



U[P] = U C [P] 



Pfr) -P(r1 



X(r) 



-dr. 



(5) 



2 When there is no possibility of confusion, we do not specify the variable for the operator V; otherwise, we 
indicate the variable by a subscript. 
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Performing a functional variation with respect to the polarization vector P, we arrive at an 
integro-differential equation defining the equilibrium polarization 

P(r) , Vf />,(r-)-V : P(rO dr , = 



X(r) 
which implies 

P(r) = X (r) J [p f (r') - V ■ P(r')] 1^4^' = X (r)E(r) . (7) 

Thus the constitutive relation for a linear dielectric is obtained as a result of functional minimiza- 
tion, with the expansion coefficient x(r) turning out to be the dielectric susceptibility. Inserting 
the equilibrium polarization ([7]) in ([5]) results in the well known expression for the total energy 
of the system: 

U = \f Mr) |73^T [MO - V ■ P(r')] drdv' . (8) 

Keeping two (or more) terms in the series (HI) introduces nonlinearity into the problem. The 
energy functional in this case is given by 

1 fPW'Pfr) 1 f \P(r) ■ P(r)l 2 
Performing a functional variation as above we now obtain 

P(r) = X (r)E(r) - 2^ [P(r) • P(r)] P(r) . (10) 

Given that the first term on the right hand side is the dominant one, we can obtain the solution 
via iteration. The first approximation would be the same as the result for the linear dielectrics. 
Substituting it back into (fTUj) . we obtain at the second order of approximation, 

P(r) = X (r)E(r) - 2^ [E(r) ■ E(r)] E(r) . (11) 

/i(r) 

One can continue with this to obtain a series of terms with higher and higher powers of [E ■ E] . 
This gives the desired result for nonlinear dielectrics. We should mention once more that this 
solution is true for weak fields so that the higher order terms are successively weaker. To ensure 
this condition we require p(r) >> x 3 ( r ) to be true to any order of approximation. 

Let us now solve ([7j) for the case of linear dielectrics. We simplify the analysis by choosing 
the (scalar) induced density p\ = — V • P as our variable. 
Using the relation 

47T(5(r - r') , (12) 



we obtain from ([7]) 

V ■ P(r) = V X (r) • J [p f (r') - V ■ P(r')] dr' + 4vr X (r) [p f (r) - V ■ P(r)] (13) 
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which implies 

e(r) Pi (r) = -V*(r) ■ j [MO + Pi(r')] dr' - 47r X (r)p f (r) , (14) 

where e = 1 + 4:irx- Equation ( fT4l) relates p; and pf . We may rewrite this equation as 

e(r)p t (r) = p f (r) - V X (r) ■ j -^^p t (v')dv> (15) 



or 



Pt(r) = ^ - -^Vx(r) ■ / -^- Pt (r')rfr' (16) 
e(rj e(r) J |r — r'| d 

This integral equation is the most general equation for total charge density in linear dielectric 

media. Note that it is a simple scalar equation for the induced charge p i; as opposed to flj[J), a 

vector equation for the polarization P whose numerical solution also requires calculation of V- P. 

Once (ED) is solved for p t , the polarization field is straightforwardly obtained by substituting p t 

for pf — VP in (JTj). The advantages of switching to the induced charge persist even in the case 

of nonlinear dielectrics. 

For a system with uniform susceptibility, we obtain the expected screening pt(r) = so 
that Pi(r) = —(1 — i)pf(r). The second term in f[T5"j) generates induced charges due to non- 
uniformity of dielectric medium. In the case of a sharp boundary, the proper limit of this term 
gives rise to surface charges. A planar interface example is described in appendix lAl 

We may rewrite ([TBI) in the form of an operator equation 

(I + C)p t = ^, (17) 
e 

where the operators I and C are defined as 

[I h] (r) = J 6{r - r')h{r')dr' , (18) 



C W = / ■ TT-i^ h ^)dr' ■ (19) 



Vx(r) r — r 
e(r) |r — r 

We will frequently make use of the kernel of this operator defined as 

C(r,r') = ^f T i^ i . (20) 
e(r) |r — r \ 6 

Note that C is completely determined by the geometry regardless of the position of the source 
charge. 

Using the formal inversion of I + C 

[I + Of 1 = I - C + C 2 - C 3 + • ■ ■ , (21) 
one may obtain the total charge density 

p t = [I-C + C 2 -C 3 + ---]^. (22) 

If the off-diagonal part C(r, r') is small compared to the diagonal delta function, series (j2lj) 
converges quickly. 



7 



III. THREE CASE STUDIES 



In this section we apply our energy minization method to three examples for which the exact 
solutions or the equations governing the exact solutions are known. 

A. A planar interface 

Let x depend only on one spatial variable z. For z > a, % = Xii an d for z < —a, x — Xi- m 
the range —a < z < a, x is a smooth function of z. Then 

c M = #-r^ = T¥^3- < 23 > 

e{z) |r — r| d e{z) \r — r'\ 6 

Let us put a free point charge q at z = d > a so that pi (r) = q5(r — dz). The total charge density 
([22]) becomes 



Pt r = -5(r-dz)-—— / —5 (r - dz )-dr 

ei 4ire{z) J |r — r\ 6 e x 

e'(z) f z — z' e'(z') z' - z" .. q , . , .. , nA . 

+ - \\ / 1 ^ - Vtt i rr5(r — dz)—drdr H , (24) 

A<ite{z) J |r - r'| 3 47re(z') |r' - r"| 3 v ; ei v 7 

where we have used e = 1 + 47IX 

In the a — > limit, e'(z) = 5(z)(ei — €2), so 

p,(r) = ^(r - dl) - -£=^f (z) J J^fHr" - di)Ur> 

+ {l^) 2 ^J0T^ ) ^ S{l "-< dl ' dl '' + --- ■ (25) 

Note that each term from the second order on has a factor of z8(z) which is zero for any z. We 
finally obtain 

*M = ^<'-*>-4^=V W i^T (26) 



The surface charge density 12J depending on the radial vector p in the x — y plane 



q 2( gl - e 2 ) d 
= 4^(e 1 + e 2 ) |p-d*|3 ' ^ 

is then obtained by setting e{z = 0) = (ei + e<i)j2. The validity of using the average dielectric 
constant at the boundary is justified by the following argument. Let there be a surface charge 
density a at the boundary. It creates an electric field of magnitude 2na directed along the normal 
vector to the surface. Assuming that there are no free charges at the interface, the boundary 
condition requires that (E± + 2na)ei = (E^ — 27ra)e 2 , where E± is a normal component of electric 
field produced by sources other than a. Therefore, a(e± + e 2 )/2 = (e 2 — ei)E±/4:TT, in agreement 
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with setting e(z = 0) = (ei + £2)/2. In Appendix lAl we present a thorough derivation of the 
a — > limit, which arrives at the same conclusion without invoking ^-functions. It is worthwhile 
to point out here that the surface charge density arises entirely from the term containing the 
gradient of the susceptibility. Our formulation is straightforward in this respect when contrasted 
with methods that first neglect the gradient of x an d then introduce a surface charge density by 
hand Q. 



B. A point charge outside of a sphere 



Consider a ball of radius a\ centered at the origin and a point charge q located at point L, 
Pf(r) = q5(r — L). In this subsection, we first obtain a set of equations for the general case of 
spatially varying susceptibility, assuming only that it changes in the radial direction. We then 
consider the case of a sharp boundary and show that the simplified expressions for the induced 
density coincide with the known results 22]. 

Let the susceptibility change in the radial direction from some value xi inside the ball to 
another value Xo outside. Gradient x is then directed radially, 



&x 

V *(r) = Qff = — r, 



and we find for C(r, r') 



C(r, r') 



Aire(r) 



e'(r) 



~/\3 



Aire(r 



-Br, 



Let us calculate 



dr'C(r, r 



q e (r ) 



e r 



-d r 



1 



{21 



(29) 



(30) 



e Q 4ire(r) [r — L| 

Assuming, for simplicity, that the point charge is located far enough from the ball, so that 
e'(r) ^ only where r < L (a generalization which would lift this condition is straightforward), 
we obtain the first order approximation for the induced charge density, 

Pi 



Pi (1) (r) 



where 

and the expansion 

^ oo £ 



c 



lm 



An q e'(r) Ir 



i-i 



21 + 1 e 4vre(r) L l+1 



(31) 



(32) 



EE 



An r 



21 + 1 rl+ 

1=0 m=-l > 



lYimir^Y^fy), r< = min(r 1 ,r 2 ), r> = max(n,r 2 ) (33) 



|ri - r 2 | 

was used. Note that any one of the spherical harmonics can bear the complex conjugation sign. 



9 



The next order is obtained by applying the operator (— C) to p^: 

p^(r) = [-C- Pi «] W = £ [/ dr'{-C{vy))p^{r')Y lm {r') 



YL{L). 



The angular integration in ()34p can be performed analytically using (129]) and (1331) : 



C ?r'(-C(r,r'))p£(r')y im (f') 



"'-ft /dr^— ^ 



r — r 



e' r 



47re(r) ' 

V ' I'm' 

The orthogonality relation for the spherical harmonics, 

J d?Y i r ml (?)Y lm (?) = 5 lll 5 n 

removes the sum, so we obtain 



dr I q dr'j^-J^r') \ dr>YZ m ,{r>)Y lm {r>) 



i in 

47r e'{r) 
21 + 1 4yre(r) 



oo I— 1 



PpP^Mdr' -(/ + 1) 



r ( r '\l+2 



o 



r 



1+2 



The same derivation leads us to a general recursive relation 
pi (n+i) (r) = J2pt +1 \r)Y lm (r)Yr m (L), 



Pt +l \r) 



lm 

An e'(r) 
21 + 1 47re(r) 



(n) 



r ( r') l+2 {n)l 



Therefore, using (I22p . we write the induced charge density for the general case of a sph 
a radially varying susceptibility as 

Pi(r) = J2pUr)Yi m (r)Yi*m(L), 

lm 

oo 

PiJj) = £p^(r), 



n=l 



where '(r) can be found via (1321) and 
In the limit of a sharp boundary, 



we immediately find that 



e'(r) = (e Q - e x )5(r - ai), 
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while the higher order contributions, 



ftt +1, M 



AlT 



n+1 



2/ + ly V 47re («i) 



"+ 1 ij-i 



la[~ 



(42) 



are found from (1381) using the generalized definition of the Dirac ^-function, 



h(x)5(x) = -MO). 



Finally, we sum all the contributions to obtain the total charge density: 

oo \ 

i+E(- c ) nW/ 



(43) 



Pt[T) 



n=l 



e Q e Q V47re(ai 



9_ (*o_ 



i-i 



E 

n=l 



e Q - ei 1_ 

' 2e(ai) 2/ + 1 



H m (f)y4(L) 



n-l 



(44) 



The sum in square brackets is a geometric series with common factor less than 1 for all I. 
Substituting e(ai) = (e Q + ei)/2 again, we derive 



Pt(r) = ^5(r - L) + ± (e - e 1 )6(r - a,) ^ Y^Y^L) . (45) 

e D ^ [{I + l)e + kij L'+i 



For the case in which the point charge is inside the ball, similar analysis leads to 

00 7,-1 J I 

Pt(r) = U{x — L) + — (e Q — ei )*(r - a,) £ + i -^Ur^(£), L < «i (46) 

/m Y\ L L ) € o + <e lJ Oi 



Using the addition theorem for spherical harmonics, 



Air 



Pi{r-L) = — i Y,Y lm (r)Y l * m {L).„ 



(47) 



m=—l 



and placing the point charge on the z-axis, L = (0,0, L), one can further simplify the derived 
equations: 



p t (r) = — d(r - Lz) + 

e Q e D 47r 



'fr-oE ktt^tm < 4f 



Pt r = — 5 r-L^ 

ei ei An 



^ (/ + 1)(2/ + 1) L' 
5(r - ax) ^ — — , , ;j , 9 P;(cos i 



^ [(/ + l)e + / ei ] a'+ 2 " 



L < ai. (49) 



where 6 is the polar angle of r. These expressions provide the correct results for the surface 
charge densities which can be found in 22 
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C. Multiple charges and multiple spheres 



We now generalize to the situation of many point charges and many spheres. In this case only 
the exact equation, not the exact solution, is known 0. According to the linear superposition 
principle, the induced surface charge on each sphere may be computed by using one free charge 
at a time and then adding up the contributions. 

Let us consider iV dielectric spheres of various radii and dielectric constants immersed inside 
a dielectric medium of dielectric constant e Q . The location of sphere % is R«, its radius is aj, and 
its interior has dielectric constant ej. No two spheres are in contact with one another. There are 
K point charges qi located at gj so that the free charge density reads pf (r) = J2f=i Qi^{ r ~ gi)- 
We assume that the variation of susceptibility in the vicinity of each sphere is radial with respect 
to the center of that sphere: 



Vx(r) 



N a N 



i=l 



i=l 



47T Ti 



(50) 



Here and throughout this section we use the tilde sign to denote radius vectors centered at the 
corresponding spheres, r = Rj + f j. 
From ffToD we have 



Pt(rJ 



Pfl r J 



e r 



47re(fj 



r-r / a . / _ Pi r 
Z — Z7^Pt( r ) dr 



e r 



(51) 



where C« plays the role of C in ( TT71) . 

Concentrating on the equation associated with a particular sphere k, we decompose /?t(r) as 



Ptlrj = Pk(r) 



Pfl r J 



e r 



(52) 



where Pi(r) is the total charge density near the surface of sphere i. Since we consider nonover- 
lapping spheres, CjCj = for i ^ j. Therefore, when focusing on a spatial point near sphere 
k, the only contribution to the overall charge density is Pk{^), so pt(r) = Pfc(r) for r sufficiently 
close to sphere k. Then in vicinity of sphere k the charge density becomes 



Pk{r) 



which may be expressed symbolically as 

[I + C fc ] p k 



r.'\3 



Pf(r 



e r 



Pk(r 



(53) 




with 



C k (r,r' 



4ire(r k 



f'|3 



(54) 



(55) 
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This implies a symbolic solution for p k 



Pk = -[l-C k + Cl-Cl + -..]C k [^ + J2 Pi ) • (56) 



Notice that the solution for the series acting on the free charges part will be essentially the 
same as that for the one sphere problem dealt with in the previous subsection. Let us consider 
Cfc Pjjtk- 

n e'(h) ~ [ r-r' , 

Ctft = ^'iFrf ft(r)dr (5?) 
We switch to vectors centered on the corresponding spheres so that the final expression is in 
terms of the local polar angle of r k , which allows easier manipulation later. In this notation, 



e'(h) 



47re(r fe ) k J \r k - (r$ - L^ fc )l " 3 

where = Rfc ~ Rj = — ^k->j represents the vector pointing from the center of sphere j to 

that of sphere k. Using the expansion (1331 . we obtain 

lm i J J i 



The angular integral in the above equation was solved by Yu [23j and employed in [9| where 
Pj oc S(fj — dj). The process for calculating pj is not affected by the detailed result of the 
integration. For now, it is sufficient to point out that the integral gives rise to a geometrical factor 
with some factorials multiplied by the multipole moment Q\ m of the surface charge distribution 
of sphere j. Denoting the integral by Aj m (dj, Lj^ fc ), 

/Y* f r j~ L j^fc \ 

we may then write 

<W') = -^E^^^'^)^^^ • (61) 

lm 

For the case of sharp boundaries between the spheres and the external medium, one then obtains 
Pi {r) = -j^AKh - a*) E aiTl [ ,fl *" lA ^' L ^ fc) ] Ylmih) ■ (62) 

lm 
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Applying the operator once again and performing the integration in the radial direction, 
we find 



C 2 n(r)- ( e -- € A 2 Kh-a k ) f 

C ^' (r) --U^K)J 2 J 



\2-2r k -r h \W£U + 



(63) 



After performing the angular integration, pj (r) becomes 



/ Co — Cfc 



47T 



.2c( flfc )(2/+i) ; 2/+T P a fc -lA im(%> L i— *)] yah) 

(64) 



It is easy to see that this process continues and one ends up having 



( € o — e fc \ 



/ e o _ Cfe 



n-l 



An 



(65) 



and therefore 



71=1 



\4ne(a k )J 



S(r k - a k ) x 



E 



lm 

e Q - efc 

47T 



E(-i)" 

n=l ^ 



/ e — £fc 



E 

lm 



2e(a k )(2l + l) 

S(r k - a fe ) x 
(2/ + 1) 1 4vr 



n-l 



47T 
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I H~ lA L( a ^h^k)]Y lm (r k ) 



_(l + l)e + le k 
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I IK lA L( a j^k)} Y lm {r k ) 



(66) 



where e(a k ) = (e Q + e k )/2 is used. We are now in a position to write down the full solution using 
(USD, (USD, and dSSI). Defining J fe = {z |a fe > |& - R fc | } and O k = {i \a k < |& - R fc | } to be the 
sets of charges inside and outside sphere k, respectively, we find 



j-i 



0* 



g» - Rfc 

- [(Z + l)e G + le k ] | gi - R fc |'+i J,m V|& - R fc I 



lm 



lm{ r k) 
Ylm(f k ) 



■E 



e Q - e k 
4tt 



5(r k - a fc ) ^2 



lm 



{21 + 1) 



_(/ + l)e G + /e A 



47T 



2/ 



j-pai-^Co^L^)]^*) (67) 



which, with appropriate rotations and taking a single point charge at the center of each sphere, 
is equivalent to (11) in fy. 
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Figure 1: Radial dependence of (a) the dielectric constant e{r/a) and (b) e'(r/a)/(47re(r/a)) for a 
monotonic step (red broken line, Eq. (|68p ) and for a non-monotonic step simulating a hydration layer 
(blue solid line). The dielectric constant changes smoothly from e± = 4 inside the sphere to e = 80 
outside. The effective radii ro = 1.13a and ro = 1.17a, respectively, are chosen so that the Born solvation 
energy in each case is equal to that in the case of a sharp boundary at radius a (shown with dotted 
line). The half- width of the steps Sr = 0.2a. 



IV. NUMERICAL CASE STUDY 



In this section we present results of numerical computations comparing the force between 
two charged identical spheres with sharp boundaries to the force between two charged identical 
spheres with smeared boundaries. For brevity, the spheres with smeared boundaries will be 
called "fuzzy spheres" and the spheres with sharp boundaries will be called "rigid spheres" . The 
dielectric constant e\ — 4 inside the spheres and e Q = 80 outside. For the fuzzy spheres there is 
an interface region ro — 5r < r < ro + Sr in which the dielectric constant changes smoothly from 
e\ to e Q in the radial direction (with respect to the center of the corresponding sphere). 

The simplest polynomial smoothly connecting e\ and e Q , i.e., satisfying the conditions e(r — 
Sr) = e±, e(r + Sr) = e D , e'(r — Sr) = e'(r + Sr) = 0, is cubic, so that the dielectric constant 
can be denned around each sphere as 



e(r) = ei, r < r — Sr 
(r-r ) 3 o r-r 



e r 



e r 



Sr 3 Sr 
e Q , r > r + Sr. 



<u - e Q ei + e Q 
4 + 2 ' 



r o — Sr < r < r + Sr 



(68) 



With a fifth order polynomial one can request additionally that e(r +Sr}i) = en and e'(r +5rH) = 
0. Letting en = 70 and Sr^ = 0.5Sr yields a non-monotonic profile, which may be used to simulate 
the hydration layer phenomenon in bio-macromolecules and clusters (see Fig. [T|). 

Let there be point charges q\ and q2 at the centers of spheres 1 and 2, respectively. The 
induced charge density is found for rigid spheres as the self-consistent solution of ( IBTl) for pi{f\) 
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and p 2 (r 2 ). Of course, (|6"7|) simplifies dramatically in the case of two spheres and two free charges. 
For fuzzy spheres, one has to use a continuous version of (|67j) in which summation over n in f[6"6l 
is carried out numerically with the n th -order terms (Itjol) calculated recursively via numerical 
integration, analogously to the method for a point charge outside a sphere, see (j38j) . ()42p and 
OH]). Notice that the / = components of the induced densities can only be produced by the 
free charge inside the corresponding sphere. Notice also that the free charges in the centers of 
the spheres induce only I = 0, i.e. spherically-symmetric, components. For these reasons it is 
convenient to distinguish the / = and I 7^ components of the induced charge density. 

In accordance with (JSJ), the total energy of the system consists of the following terms: 

(i) interaction between the point charges (screened by ei), 

||, (69) 

where L is the length of the vector Li_> 2 — — L 2 _ +1 , connecting the centers of the two spheres, 

(iia) interaction between each point charge and the I = component of the induced charge in 
the interface region of the other sphere, 

(iib) interaction between each point charge and the / 7^ components of the induced charge 
in the interface region of the other sphere, 



2 V. J |r 2 + Li_ 2 | J |ri + L 2 -i| 



(71) 



(iiia) interaction between each point charge and the / = component of the induced charge 
in the interface region of the same sphere, 

^ ?l |M^ ifl + ?2 |M|U if2 ^ (72) 

(iiib) interaction between each point charge and the / 7^ components of the induced charge 
in the interface region of the same sphere, 

The sum of terms (i) and (iia) is equal to the energy of interaction of two point charges in 
dielectric medium e Q 

2*. (74) 
e L 

This energy is the same for rigid and fuzzy spheres. In contrast, terms (iib) are different for rigid 
and fuzzy spheres and are the main source of differences in the forces in these two situations. 
Finally, terms (iiib) are zero for the point charges located at the centers of the spheres, while 
terms (iiia) are the Born solvation energy in this case. 
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Figure 2: Radial dependence of the induced electric density p(r/a)\ l=0 for the monotonic (red broken 
line) and non-monotonic (blue solid line) steps shown in Fig. [TJ The density is normalized by the value 
of the point charge in the center of the sphere. The inset magnifies a small, oscillatory feature associated 
with the non-monotonic step. 

Born solvation energies are quite different for rigid and fuzzy spheres, since for fuzzy spheres 
the induced charge density tends to accumulate near the inner boundary of the interface region. 
Indeed, the operator C is proportional to e'(r)/e(r) and e(r — Sr) = t\ <C e(r + Sr) = e Q . 
This asymmetry is present at each order n and is preserved after the summation over n. Radial 
dependences of the I = components of the induced densities are illustrated in Fig. [2j On 
the other hand, fuzzy and rigid spheres model the same physical objects, so it is reasonable 
to assume that whatever profile of the dielectric constant is chosen, the Born solvation energy 
should remain the same. For this reason, we adjust the effective radius r for each profile of the 
dielectric constant so that the Born solvation energy is equal to that of a rigid sphere of unit 
radius, see Fig. CD 

In Fig. [3] we present the dependence of the interaction energy on distance for a pair of rigid 
spheres and for two pairs of fuzzy spheres, with monotonic and non-monotonic behaviour of the 
dielectric function in the interface region, respectively. The energies are normalized to the energy 
of interaction of point charges (1741) . The forces between two fuzzy spheres and between two rigid 
spheres are shown in Fig. HI The forces are normalized by the interaction force between two point 
charges. We note that the seemingly weaker effect for the fuzzy spheres with non-monotonic e(r) 
dependence is due to the fact that e(r) changes faster near the inner surface of the interface 
region to make room for the feature representing the hydration layer. This makes the fuzzy 
spheres with non-monotonic e(r) dependence effectively more similar to rigid spheres for fixed 
Sr (compare the charge density distributions in Fig. [2]). 

For very thin interface regions (Sr — > 0), the forces between two rigid and two fuzzy spheres 
are equal, as expected. For fuzzy spheres with moderate interface region widths, the repulsion 
increases with the width. However, this trend quickly saturates (Fig. [5]). Qualitatively, this 
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J I i I i I i I i L 

2 3 4 5 6 7 
Ua 



Figure 3: Energy of interaction between two spheres with sharp (thin line) and smeared (thick lines) 
boundaries. The red broken thick line corresponds to the case of the monotonic radial dependence of the 
dielectric constant, while the blue solid thick line corresponds to the non-monotonic radial dependence 
shown in Fig. [TJ Free charges of the same sign are located at the centers of the spheres. The energies 
are normalized by the Coulomb energy of these point charges in the uniform dielectric medium e . The 
vertical dotted lines indicate the contact points. 




u i n i , i , I i I i i 

2 3 4 5 6 7 
Ua 

Figure 4: Interaction forces between two spheres with sharp and smeared boundaries. The line identi- 
fications are same as in Fig. [3l 

saturation can be explained by two opposing effects. The increase in the interface width increases 
the size of the spheres thereby strengthening the repulsion. On the other hand, the induced charge 
density tends to concentrate near the inner surface of the interface which remains around r = a 
to maintain constant Born solvation energy. Therefore, the bulk of the induced charge on one 
sphere becomes farther from that of the other sphere, hence weakening the repulsion. 

We finally note, that if the point charges are located away from the centers of the spheres, the 
terms (iiib) depend on the relative position and orientation of the spheres. In this case one can 
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Figure 5: Maximum difference in interaction forces between two spheres with smeared and two spheres 
with sharp boundaries, occuring at the contact point 2(ro + Sr), as a function of half width of the 
interface region 5r. The forces are normalized by the interaction force between two spheres with sharp 
boundaries. The red broken line corresponds to the case of the monotonic radial dependence of the 
dielectric constant, while the blue solid line corresponds to the non-monotonic radial dependence shown 
in Fig. [U 

still define the Born solvation energies as the sum of terms (iiia) and (iiib) at large separations, 
but the terms (iiib) would contribute to the difference of interaction forces/energies between the 
rigid and fuzzy spheres. 



V. CONCLUSIONS 



We have presented an energy minimization formulation of electrostatics that allows compu- 
tation of the electrostatic energy and forces to any desired accuracy in a system with arbitrary 
dielectric properties. We have derived an integral equation for the scalar charge density from 
an energy functional of the polarization vector field. This energy functional represents the true 
energy of the system even in non-equilibrium states. Arbitrary accuracy is achieved by solving 
the integral equation for the charge density via a series expansion in terms of the equation's ker- 
nel, which depends only on the geometry of the dielectrics. The streamlined formalism operates 
with volume charge distributions only, not resorting to introducing surface charges by hand as 
is done in various other studies of electrostatics via energy minimization. Therefore, it can be 
applied to arbitrary spatial variations of the dielectric susceptibility. The simplicity of applica- 
tion of the formalism to real problems has been shown with three analytic examples and with a 
numerical case study. We found that finite boundary widths introduce a measurable correction 
to the interaction forces as compared to sharp boundary case. For two charged identical spheres 
the correction is about 10%. 

The formalism has various potential applications in modeling electrostatic interactions be- 
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tween solvated molecules: it enables one to go beyond the widely used simplification of atoms 
and molecules as dielectric balls immersed in a dielectric solvent, as was first suggested by Born in 
the early twenties of the last century 24|. For example, the description of an aqueous solvent as a 
continuous and homogeneous dielectric medium fails to account for the strong dielectric response 
of water molecules around charges. Normally, charged ions and surfaces give rise to hydration 
layers by orienting and displacing surrounding water molecules. These hydration phenomena 
are very important in many biological processes such as protein folding, protein crystallization, 
and interactions between charged biopolymers inside the cell. With our formalism one can now 
consider arbitrary structures for such hydration layers and arrive at a possibly more realistic and 
reliable analysis of the molecular mechanisms in bio-chemical interactions. 

Applied to MD simulations, this formulation is still an implicit solvent scheme, and the 
position-dependent susceptibility is therefore a model parameter (indeed, the only one). To 
obtain an estimate of the macroscopic dielectric susceptibility at the molecular level or at the 
intermolecular boundaries one has to explore physics at the atomic level and introduce some 
coarse graining. Given that the dielectric susceptibility is related to the charge fluctuations 
as a response to external perturbations, one can estimate susceptibilites through the study of 
linear/nonlinear response. For example, the dielectric susceptibility can be related to the cor- 
relations of the net system dipole moment and local polarization density 25j. A fully quantum 
mechanical treatment of solvation of biological systems might be hindered by limits of numerical 



accuracy [26] and will demand much more computational power than currently available. We 
believe that quantum mechanics, in particular, density functional theory, can in principle be 
used to calculate the local dielectric susceptibility which in turn should be used as input for the 
implicit solvent methods, such as the one described in this paper. 
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Appendix A: SHARP BOUNDARY LIMIT IN THE PLANAR INTERFACE PROB- 
LEM 

Let us demonstrate how a rigorous limiting procedure applied to fl25l) produces correct ex- 
pression for the surface charge density in the case of sharp planar interface. The surface charge 
is found by integrating the charge density over the range —a < z < a in which \ changes from 
X2 to xii an d then taking the limit a — > 0. 

We return to (12^1) and, making use of the azimuthal symmetry of the problem, expand the 
kernels in terms of Bessel functions J m |l2f , 
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1 



|r — r 
1 



£ 

m=—oo ' 

"CO 



J m (kp) J m {kp')e- k ^>- z <Mk 



d z\ 



J (kp)e- k{d - z) dk. 



(Al) 
(A2) 



Here the position vectors r and r' are represented via the polar vectors p and p' in the z = 
plane, r = p + zz and r' = p' + z'i. The polar vectors are in turn defined through their lengths 
p = \/x 2 + y 2 and p' = a/ a;' 2 + y' 2 and their azimuthal angles and 0'. The notation z > (z<) is 
used for the greater (lesser) of the corresponding z and z'. 

We now treat each of the terms in the expansion of separately. The first term is the 
screened point charge. All other terms form the induced charge density at the interfacial region. 
The first contribution to the induced charge density is given by 



Pi (1) (r) 



q e [z) 



z — d 



e\ ^Ke{z) |r — dz\ 3 



The corresponding surface charge density is 



^ (1) (P) 



-— lim 



a e'(z) 



-d 



\p — dz\ 3 



+ 0(z) 



dz. 



All the O(z) terms vanish since for any bounded function h(z) 

/a pa 
z n h(z)dz < lima™ / \h{z)\dz = 0, V n > 0. 
J-a 



Thus, 



^ (1) (P) 



q 



d 



(A3) 



(A4) 



(A5) 



(A6) 



4nei \p — dz\ 

Here we have used the notations f(z) = hi[e(z)], fi = f(a) = ln[ei] and f 2 = f{—o) = ln[e2]. 

We can similarly evaluate all the other contributions to the induced surface charge density. 
The second contribution to the induced charge density is 



Pi (2) (r) 



z-z' e'(z') 



z'-d 



€i 4-Ke(z) J |r — r'| 3 4ire(z') \r' — dz\ 



: p'dp'd(f)'dz' 



Using (]Aip . flA2j) . and the completeness relation for Bessel functions 12| . 

f°° 1 
/ J m (kp)J m (k'p)pdp = -5(k - k), 
Jo k 

we obtain, after integration over <$' and p', 



(A7) 



(A8) 
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Pi (2) (r) 



q e'(z) d f e'(z') 



ei 4:7ie(z) dz J Ane^z' j j 



J (kp)e~ k{d ~ z ' ) e~ k(z> ~ z<) 2iidkdz' 



a e'(z') 



eu 



dz' 



2 e'(z') 



-2k {z - z >) dz , 



(A9) 



The corresponding surface charge density is then 



./'u'tf^! r 

ei^oJ_ a A7re(z)2J 



kdke- k{d - z) J (kp) x 



a e'(z') 



eu 



2 e'fz'l 



-2fc(z-^) dz / 



Applying to (lAlOj) the same argument used in deriving (1A6|) . 



0"i 



.(2) _ I 



/ kdke~ kd J (kp) lim / dz- 
Jo J -a ' 



e'(z) 1 



ei Jo a ->"J~a 4vre(2)2 

The integral over k is evaluated using ( 1A2I) as 

d 



a e'fz'] 



eu 



dz' 



2 e'(^ 



eu 



dz' 



U (kp)e- kd dk = I J (kp)e~ k{d ~ z) dk 



Then 



.(2) 



q d 



lim / dz 



z=0 



e U 



dz |r — d z| 



d 



2=0 



|p — dz| 



ei |p — dz| 3 a^o J_ a Aixe(z 
Finally, we obtain that a^ 2 ' = 0, 

d rf 



gC/i + /a) - /(*) 



0"i 



.(2) 



_J? 

47rei |p — dz| 



df 



(fi + / 2 ) - f{z) 



0. 



(A10) 



(All) 



(A12) 



(A13) 



(A14) 



Analogously, the expressions for the induced surface charge densities up to the fifth order are 
found to be 



(i) 



(2) 



(3) 



(1) 



(5) 



q d 



;(/l-/ 2 ), 



Atxei |p — dz| 3 
0, 

1 d rV3 



4ne 1 |p - dz| 3 12 



0. 



q d 



Ane 1 |p- dz| 3 120 



(fl - / 2 ) J 



(A15) 
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In general, the surface charge density is of the form 

„<»>(;) = -'ton rdz /W Z ~ d -. 

e 1 a^o J _ a Ane{z) |r — dz\ 6 

1 



a e'(z' 



eiz 



z e'(z" 



'^ n - l \f{z'))dz' 



ez 



- — lim - 



dz 



e'iz) z — d 



ei a^o 2 J_ a 4ire(z) |r — dz\ 



Q 



d 



g {n \f)df. 

Aixei \p - dz\ d J f2 



fi 



The functions g^ n '(f{z)) up to n = 5 are 



(A16) 
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^(/i + /a) 




{2) U\z)) = 


-m -+ 




{3) (/w) = 


2 


^(/i + / 2 )/(*) 




(4) (/w) = 








6 


{5) (/w) = 


/ 4 (^) 
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^(/l + / 2 )/ 3 ( 


^) + ^i/ 2 / 2 W + ^(/i + / 2 )(/ 1 2 



-4/i/ a + #)/(*) 



■^/i/ 2 (/i -3/i/ 2 + / 2 2 ). 



(A17) 



We will show by induction that g^ n \f) is 

s (n) (/) = (- 



1 1 

(n-l)! J 2 



Cxg^if) - ^C 2 g^ 2 \f) + ic 3 ^ (n - 3) (/) + 



\n— 1 fn—1 ^ " 1 / -j^n— m— l£f 



_-j\n— 1 yr 



(n- 1)! 



m=l 



(rt — m)! 



- _.(m)/j\ 



(A18) 



where the coefficients C n = + f£- First, (1A18|) can be explicitly verified up to n = 5 using 
(1A17p . Second, we show that if this expression holds for some integer n, then it also holds for 
n+ 1. From Eq.( 1A16j) we can write, 
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9 {n+1) {f{z)) 



fi 



/(*) 
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m=l 



(n — m)! 



./l 

n-1 



ff (m) (/)4T 



/(*) 



/2 



(-i)-r ^ i (-ir- i (/r + / 2 n ) + 1 ^2 (- 1 ) w " m " 1 c , n- 



/>,! 



1(1 = 1 



(n — m)! 



9 {m) U)df 
(f) 



j(n+l)-l 

((n + l)-l)! 



(n+l)-l 



+ ^ E 



m=l 



/ i \(n+l)— m— 

( J » -r, {n ^~ m 9 {m \f) (A19) 
((n + 1) — my. 



We thus proved that g^ n \f) is given by (1A18[) for any given integer n > 2 with g^\f) = 1. 
We now need to find the integral J <t>) in (IA16p . We will show by induction that 

•A . . £ 

(A20) 



/2 



(») (/) = _ 2 ^. 

n! 



where w = /i — /2 and are the coefficients of the expansion 



e u + 1 ^ n 

n=0 



(A21) 



It is easy to see that Eq = 1. 

The base for the mathematical induction for ( 1A20I) is easily established for the first few terms 
using ( 1A17I) . Now we verify that ( 1A20I) holds true for n + 1 if it is true for n. To do so, we 
integrate both sides of flA19j) and use the assumption flA20j) to obtain 
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m! 



(n + 1)! 



In the second step we have included an m = term in the summation and in the third step we 
have added and subtracted an m = n + 1 term. It can be easily verified that the right hand side 
of ( 1A22I) is the s n+l term of the following expression. 



- 2e~ /lS + [e 



-fis 



-fas 
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This completes the proof. 

Summing over all the terms, we have 



-A n nl V e u + lj ei + e 2 

n=l J- 3 n=0 



We note that the series converges for |u| = lne /ei < ir. This means that if one medium is water 
(e Q ~ 80) then for the other material the dielectric constant e± > « 3.47. However, using 

techniques similar to Borel summation, one can show that the series can still be summed to the 
correct final formula for larger values of \u\. 

Finally the induced surface charge density becomes 

q 2(ei-e 2 ) d 
v ' 4vrei ei + e 2 \p-dz\ 3 K J 

which is identical to ( 1271) . Thus, we have rigorously justified using the average dielectric constant 
(ei + £2)/2 at the boundary. 



Appendix B: EVALUATION OF A FOR SPHERES WITH SHARP BOUNDARIES 

To compute A J lm (aj, L,-_>fc), defined as 

* ( r j" L j^fc 



AL(a„ L^ fc ) ee / ^^'-^ > ( Bi) 

J r j _L, j-»fc 



for the case of spheres with sharp boundaries, expand the charge density on sphere j as 



to find 



lmK 3 3 ' lm J L +1 fc (l + t 2 -2tcos^)( z + 1 V 2 V 3 31 3 ^ 1 



where use has been made of the geometrical fact that \ij — ~Lj->k\ = %-»feV 1 + ^ 2 — 2tcos9j with 
t = fj/Lj^k- The delta function renders the radial integration trivial: 



where d and are the polar variables of (fj — ~Lj_>k)/\ij — Lj— >fc I and t = aj/Lj^ now. All of 
the angular variables are measured with respect to a coordinate system whose z axis is parallel 
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to Ly_>fc. The angles $ and tp must be expressed as functions of the integration variables 9j and 



cos$ 



(tcos9j — 1) 



1 +t 2 - 2t cos Oj 

ip = <j)j . 

Since the definition of the spherical harmonics is 
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The integration over <pj produces 27r<5 mm /. The integration over cos#j is then the integral calcu- 
ated by Yu 23| . The final expression for A is 



A L( a i» L i^fc) 



QL^(-i)^ m (^ + 0!v / 2T+T 



^ tig k [4>ir(l + m)\{V + m)!(i - m)\{V - m)\(2l> + l)} 1 / 2 



(BIO) 



where Q;, m = 47ra|cr/, 
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